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Abstract 

In this paper, the nonlinear mode conversion of extraordinary waves in nonuniform magnetized 
plasmas is studied using the variational symplectic particle-in-cell simulation. The accuracy of the 
nonlinear simulation is guaranteed by the long-term accuracy and conservativeness of the symplectic 
algorithm. The spectra of the electromagnetic wave, the evolution of the wave reflectivity, the 
energy deposition profile, and the parameter-dependent properties of radio-frequency waves during 
the nonlinear mode conversion are investigated. It is illustrated that nonlinear effects significantly 
modify the physics of the radio-frequency injection in magnetized plasmas. The evolutions of 
the radio-frequency wave reflectivity and the energy deposition are observed, as well as the self¬ 
interaction of the Bernstein waves and mode excitations. Even for waves with small magnitude, 
nonlinear effects can also become important after continuous wave injections, which are common 
in the realistic radio-frequency wave heating and current drive experiments. 
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I. INTRODUCTION 


Radio-frequency (rf) wave injection plays an important ro 


e in auxiliary heating and 
4], A typical method for rf 


current drive for modern magnetic confinement fusion devices jl- 
wave injection into overdense plasmas, where the plasma frequency UJ pe is much larger than 
the electron cyclotron frequency oo ce , is to excite the electron Bernstein waves (EBWs). Two 
different approaches are commonly used for coupling rf power to EBWs, i.e., the ordinary- 
extraordinary-Bernstein (O-X-B) and extraordinary-Bernstein (X-B) mode conversion. The 
O-X-B scheme is usually used for the electron Bernstein wave current drive (EBCD) in 


spherical tokamaks (S 
method is ineffective 


V), where the electron cyclotron resonance current drive (ECCD) 
5]. The X-B scheme requires steeper plasma density profile, which 


is also available in some STs j§-9|. During the X-B mode conversion, the electromagnetic 
wave is initially injected into the plasma perpendicularly to the background magnetic field, 
and propagates as an extraordinary (X) wave within the plasma. When the X wave reaches 
the cutoff location, where the local right-handed cutoff frequency c Or is equal to the frequency 
of the incident wave uj, part of the energy of the incident wave penetrates the cutoff layer 
and continues propagating as EBWs. The X-B conversion will be significantly enhanced if 
the electron density gradient becomes large enough. 

The linear X-B mode conversion can be theoretically investigated using the cold plasma 
model & 11]. In this model, the kinetic effects of the electrons are neglected. At the same 
time, the ions offer a positive charged background. The brief description of the linear theory 
is presented in Appendix [A] The linear theory predicts that the reflectivity of the injected 
wave is a constant that does not change with time. 

The cold plasma assumption is only accurate when the wave length is much larger than 
the electron gyroradius. As the electron temperature increases, the kinetic effect becomes 
significant and should be taken into account. The kinetic theory of mode conversion is 
so complicated that the analytical solution is unaccessible except for some special cases 
js]. In addition, the linear condition no longer holds for large-amplitude or long-term wave 
evolution. In modern magnetic confinement devices, the power of the rf wave injection is 


so large that the nonlinear effects can not be neglected. The nonlinear effects durin 
propagation of Bernstein waves have been directly observed in some experiments [4, 


lie 


12 


13]. Those nonlinear effects can modify the reflectivity during the wave injection, which 
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cannot be treated by the linear theory. The simulation study on the nonlinear physics of rf 


waves in plasmas has become an active research held 


14 


17]. However, the analytical and 


numerical studies on nonlinear rf physics are severely constrained by existing mathematical 
and algorithmic tools. 

Containing both kinetic and nonlinear effects, first-principle particle simulation is re¬ 
garded as a necessary and powerful tool for dealing with the nonlinear rf physics. Different 


aspects o 


codes 
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the X-B mode conversion have lately been investigated using particle-in-ccll (PIC) 


17]. For nonlinear mode conversions, where the timescale of the mode conversion 


is much larger than the periods of the rf waves, long-term simulations are inevitably required. 
Since the number of time steps in a long-term simulation is large, the coherent accumulation 
of the numerical error from each time step turns out to be the most serious difficulty. To 
ensure the reliability of the simulation results, long-term numerical accuracy and stability 
should be guaranteed. To overcome the difficulty, a newly developed variational symplectic 
PIC method jlij is employed in this paper to study the nonlinear X-B mode conversion. The 
symplectic PIC algorithm bounds the global numerical error for an arbitrary long time by 
preserving the discrete symplectic structure of the original Lagrangian system. Its long-term 
accuracy, stability, and conservativeness have been theoretically discussed and practically 


verified 


18 


12 ]. 


Suffering no accumulation of the global errors, the symplectic PIC method shows its 
unparalleled advantage in the long-term simulations and provides accurate nonlinear results. 
In the simulation of nonlinear rf mode conversion, it is illustrated that the nonlinear effects 
significantly modify the physics of the rf wave injection in magnetized plasmas. The wave 
reflectivity, which is regarded as a constant in linear theory, varies in the nonlinear mode 
conversion process. The evolutions of the rf wave reflectivity and the energy deposition 
are observed, as well as the self-interaction of the Bernstein waves and mode excitations. 
The spectra of the electromagnetic wave, the energy deposition profile, and the parameter- 
dependent properties of radio-frequency waves during the nonlinear mode conversion are 
also investigated. It is found that even for waves weak in magnitude, nonlinear effects can 
become significant after continuous wave injections. When simulating rf waves with small 
amplitudes, nonlinear phenomena become prominent after several hundred wave periods. 
New modes are excited over a sufficient long time even if the amplitude of incident wave is 
small. To study the parameter dependency of the nonlinear X-B conversion, we compare 
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the rf wave injection processes with different electron temperatures, amplitudes of injection 
waves, and plasma density gradients. It is discovered that a lower electron temperature can 
enhance the reflectivity, and the reflectivity depends on the plasma density gradient in a 
sensitive way. 

The paper is organized as follows. In Sec. [hQ we introduce the physical model and the 
problem of the nonlinear X-B mode conversion. The numerical method used in the long¬ 
term nonlinear simulation, i.e., the variational symplectic PIC algorithm, is also introduced. 
In Sec. IIII1 the symplectic PIC simulation of the nonlinear X-B conversion is carried out. 
Different factors, which effect the nonlinear rf wave processes, are analyzed in detail. Section 
m is a the brief summary and future research plan. 


II. THE MODEL OF X-B CONVERSION AND THE SYMPLECTIC PIC METHOD 

In the X-B conversion process, the incident wave is perpendicular to the background 
magnetic field. This problem can thus be described by a planar model. Consider a planar 
plasma with the density nonuniformity in the ^-direction and a uniform background magnetic 
field along the ^-direction. The wave source is located outside the plasma. The source of 
the electromagnetic wave, with the electric field component along the ^-direction, excites an 
X wave in the plasma boundary. The X wave propagates in the plasma until it reaches a 
cut-off location, and part of the wave converts to the Bernstein wave. There are two main 
restrictions in the PIC simulation of the X-B mode conversion. Firstly, the wave length of 
the Bernstein wave is much shorter than the electromagnetic wave with the same frequency 
in the vacuum. So the grid length in the simulation should be small enough to resolve 
it, which means the time step is restricted to a small value correspondingly. A large time 
step will lead to numerical instability because the particles can pass through multiple grids 
within a single time step. Another issue is that the group velocity of the Bernstein wave 
is much slower than the speed of light in the vacuum, so the mode conversion process lasts 
for a long time. It takes hundreds of periods before the wave reaches the core region of 
plasma. Then at least ten thousand time steps are required to represent a full X-B mode 
conversion process. Traditional numerical methods may fail to obtain useful results due to 
the accumulation of numerical errors. In order to study the mode conversion process with 
more accuracy, a variational symplectic method, which has excellent long-term conservation 
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properties [19j], is adopted in the particle-in-cell simulation. 

The basic idea of constructing a variational symplectic scheme is to discretize the La- 
grangian and perform the discrete variational principle 20]. The starting point is to consider 
the self-consistent system consisting of charged particles and electromagnetic fields. The La- 
grangian is the integral 

£d 3 x , 


where 


£ = 


£o 

2 L 


£(x s , A,(j>;t) = 


-A (x, t) — V0 (x, t) 



( 1 ) 


1 [V x A (x, t)] 2 


2fi 0 


N 


Y. s 


X - X., 


s =1 


\m s ± 2 s + q s A (x, t) • x s - q s cp (x, t) 


( 2 ) 


is the Lagrangian density. Here m s , q s , x s , and ± s denote the mass, electric charge, velocity, 
and position of the s-th particle respectively, eo and /io are the permittivity and permeability 
of free space, and A and 0 are vector potential and scalar potential of the electromagnetic 
fields. The time derivative of vector potential is denoted as A. 

We follow the method used in Ref. 18] to discretize the Lagrangian L presented in Eq. (JT]) 


in cubic meshes and perform the discrete variational principle to obtain the numeric scheme. 
The discrete Lagrangian is 

Ld(l, l + 1) — Lf + Lj + L, p , 


where 


L f = 


Lp 


Lr = 



\ I eo ( A I x, (/ + £)A t 


(3) 


— fv x A fx, (/ + i)At^ ) dx > ( 4 ) 


ho 




n 


~ x > a ) m ^ x ^ dx = 5Z ’ 

S S 

y q s S(^ s ,i ~ x, a) (A (x, l At) • x Sii ) dx , 


(5) 

( 6 ) 


A | x, (7 + )At j 'y ^ (AT A fl (x Xi Jifc ) , 


i,j,k 


A (x, /At) = y Aij,k,iW (x - Xjj, fe ) , 


i,j,k 


A ( X, (/ + ^)At j = y (Aij ifci i + 1 - Aj J, fc ,i) W (x - Xj Jjfc )/At 
/ i,j,k 

± St i = (x^+1 - x SjZ )/A t , 


(7) 

( 8 ) 

(9) 

( 10 ) 
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( 11 ) 


W(x) = 


(1 — |x|/Ax)(l — \y\/Ax)(l — |z|/Ax) , |x|, \y\, \z\ < Ax , 

0 , otherwise . 

Here Ax and A t are the grid size and time step respectively. The subscript l stands for 
variables in the Zth time interval, and i,j,k are integer indices for three directions of the 
Cartesian coordinate frame. The electric potential 0 is set to be zero because of the adoption 
of temporal gauge. The corresponding discrete action integral A d is 


Ad — AtL d (l, l + 1) , 


( 12 ) 


and the motion equations are obtained from the discrete variational principle 

dA d 


dA d 


<9x 


= 0 , 
= 0 . 


(13) 

(14) 


S,l 


The discretization of Lagrangian and the discrete variational principle lead to a PIC iter¬ 
ation scheme. Theoretically, the variational symplectic PIC scheme automatically conserves 
a discrete symplectic two-form 


20 


21[. For any discrete Lagrangian L d (qi, qi+i), where qi 


denotes the generalized coordinates discretized at the Zth time step, there exists an discrete 
symplectic two-form defined as 

dL d (q h qi + 1 ) 


fly-i-i — 


-d qi A q i+ i , 


(15) 


dqidqi+i 

where the operator d denotes the exterior derivative. The discrete Lagrangian of Vlasov- 
Maxwell system is defined in Eq. (??) with qi = x Si {). The two-form in Eq. 

is a geometric structure, called the discrete symplectic structure. Because the original self- 
consistent Vlasov-Maxwell system is Hamiltonian, it is naturally equipped with a continuous 
symplectic structure, which is conserved exactly during the time evolution. The conserved 
discrete symplectic structure in the variational symplectic PIC method is an approxima¬ 
tion of the original continuous geometric structure. The good long-term properties of the 
symplectic algorithm are guaranteed through the conservation of the discrete symplectic 
two-form, that is hh-yz = for any l, during time advance in the numerical simulation, 

which can be proved exactly 20]. Therefore, the numerical method has long-term conser¬ 
vativeness and accuracy, which are necessary when simulating long-term nonlinear rf wave 
problems and have been verified numerically Q. 
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To decrease the numerical noise caused by the lack of sampling points, the shape function 
is applied. In this paper it is chosen as 


S(x;a) = Si(x; a)Si(y; a)Si(z; a) , 


where 


5i(:c; a) = aq 


fl. 

M 

1°' 

M 


x\ < a , 
x\ > a , 


(16) 


(17) 


and ot\ = l/(2a). The factor a is set to be 0.25Ax. 

With the variational symplectic PIC method, we are able to simulate the nonlinear process 
of X-B mode conversion. The plasma is inhomogeneous in x direction, and as an example 
we assume the electron density profile is 

( exp (— [x/Ax — (n r + 320)) 2 /(2(n r /3.5) 2 )) , 0 < x < (n r + 320) Ax , 

n e (x) = n 0 < (18) 

I 1 , (n r + 320), Ax < x < 1300Aa; , 

where n r Ax is the thickness of the plasma boundary, Ax = 2.773 x 10 -5 m, and n 0 = 2.0 x 
10 19 rri -3 . The density prohle is plotted in Fig. [U The physical domain is 0 < x < 1300Ax. 
At x = 0 we use the Mur’s absorbing boundary condition (MABC) 22] to eliminate the 
unwanted wave reflections. However at x = 1300Ax the density of plasma is not zero and 
the MABC fails. To construct an absorbing boundary, we apply an auxiliary calculation 
domain 1300Ax < x < 1792Ax with the decreasing electron density 


n e (x) = no exp (-(r/Ax — 1300) 2 /26122) 


(19) 


During the simulation, the Bernstein wave propagates in the plasma but not fast enough 
to reach the boundary position x = 1300A(c. At x = 1792Aa: the MABC is adopted. The 
full simulation domain is the 1792 x 2 x 2 cubic mesh, and periodic boundary conditions 


-2 


m 


are adopted in both y and ^-directions. The simulation domain is small (about 4 x 10 
compared with the fusion device, so the external static magnetic held can be treated as a 
constant. We set Bo = -BoW with Bq = 0.55T. The size of the plasma boundary and the 
electron temperature are set according to experimental data [3, 24] in National Spherical 
Torus Experiment (NSTX), where the n e /Vn e is about 0.01 m and the electron temperature 
T e is about 100 eV. We set n r in Eq. (fT8l) to be 380 in all examples except for cases stated 
otherwise. The electron temperature T e is chosen to be different constant values independent 
of spatial positions in different problems, and the distribution of electron velocity is set to 
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FIG. 1: Density profile in X-B conversion simulation examples. The Ur and ujl are the 
right and left cutoff frequencies for X wave. Here n r = 380 and to = 1.62cu ce . 


be Maxwellian. The wave source is located at x = 4Ax. It excites a linearly polarized 
electromagnetic field as 


E(4Ax, t) = E\e y sin (c ot) , 


( 20 ) 


where E] is the amplitude of the oscillating electric field, and u is the frequency of incident 
wave. The time step is set to At = Ax/c ~ 9.24 x 10 -14 s. According parameters, the 
corresponding electron cyclotron frequency is o; ce = 9.67 x 10 10 rad/s, and the plasma 
frequency is oo peo = 25.2 x 10 10 rad/s in the core plasma. The frequency of the incidence 
wave uj is chosen to be cu = 0.0145/At 1.62cu ce . According to the linear mode conversion 

theory for extraordinary waves (see Appendix [A]), the wave absorption efficiency is high 
with the above parameters. To suppress the numerical noise in the full-f PIC simulation, 
the number of sampling points should be as large as possible. However, the number of 
sampling points is limited by the computing power available. To satisfy the constraints of 
both the numerical accuracy and the computing power, the number of sampling points per 
grid at the core plasma is set to be 4000 in the numerical cases. With these methods and 
parameters, the nonlinear X-B process is numerically studied the next section. 
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FIG. 2: The contour of E y (x,t ) in the time interval [0,115/cu] 

III. SIMULATIONS OF X-B MODE CONVERSION 

A. Full wave solution in nonlinear X-B mode conversion 

The nonlinear X-B mode conversion can be clearly demonstrated by the evolution of 
E y (x,t). The amplitude of the wave source E\ is set to 1 MV/m, and the temperature of 
electrons is set to 57.6 eV. The evolution of E y in different time intervals are recorded as 
contour plots in Fig. [2] and Fig. El where the horizontal axis is the location x and the vertical 
axis is the time t. At the beginning of the wave injection, the mode conversion structure 
forms in a short time, and the reflection is stable (Fig. [2]). After about 100 periods, the wave 
structure to the right of UHR becomes unstable, as marked by the dashed ellipse in Fig. |3l 
The wave structure varies with time irregularly, which is an implication of the nonlinear 
effect. There is no direct evidence of nonlinear generation of other frequencies in Fig. 3. 
Further results of nonlinear generation of EBWs with other frequencies will be discussed in 
Sec. UlFBl 

B. New mode generation during X-B mode conversion 

The primary signature of nonlinear processes is the generation of additional modes with 
frequencies other than the incident frequency. In this subsection, the new mode generation 
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FIG. 3: The contour of E y (x,t ) in the time interval [530/a;, 645/a;] 


during the X-B mode conversion is investigated. The total number of simulation time steps 
is set to be 84000. During the mode conversion, the evolution of spectrum of E x is recorded 

in Fig. [H 


In the initial stage of wave injection, see Fig. 4(a), there only exists the wave with 


frequency u. After about 116/a;, see Fig. 4(b), the electrostatic wave propagates to the 


position at x ~ 608Ax, where a mode with the frequency 2u appears. In Figs. 4(c) and 


4(d), the component of the second harmonic mode becomes stronger, meanwhile another 

we can see 


mode with the frequency a/ = 0.77a; appears. Comparing Figs. 4(e) and 4(f) 


that the mode with frequency a/ splits into two packets, propagating inwards and outwards 
seperately. 


The generation of the mode at frequency 2a; can be explained as the self-interaction of 


the initial Bernstein wave. This second 

been studied extensively 13, H, Q 


rarmonic generation (SHG) of Bernstein wave has 


26]. Then the wave with frequency a/ = 0.77a; is 


generated near the location where the SHG of Bernstein wave happens, which implies that 
the generation of the wave with frequency a/ is closely related to the SHG of Bernstein wave. 
According to the theory of parametric decay instability (PDI), there should be a nother wave 
with frequency co", which satisfies either 2a; + a/ = to" or 2a; — a/ = co". To verify this PDI 
process, we plot in Fig. [5] the dispersion relation of electromagnetic EBWs in hot plasma, 
where the electron density of the plasma is n e (608Ax), and the electron temperature is 
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(a) [0,174/a;] 


(b) [116/a;, 290/a;] 
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(c) [290/a;, 464/a;] 


(d) [406/a;, 580/a;] 




(e) [580/a;, 754/a;] 


(f) [812/a;, 986/a;] 


FIG. 4: The evolution of frequency spectra of E x (x,t) during the nonlinear X-B mode 
conversion. In the simulation Ei = 1.0 MV/m, n r = 380, and the spatial domain of the 
simulation lies in the interval [240Ax, 960Ax]. Each sub-figure denotes the spectrum at 
different time intervals. We can see that after about 116/a; the wave with the frequency 2a; 
is excited, and after 290/a; another wave with frequency 0.77a; appears. 
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FIG. 5: The dispersion relation of electromagnetic EBWs (blue dots) and wave vectors of 
the coupled modes (solid arrows). The plasma parameters are the same as those we used 
in the X-B mode conversion simulation at x — 608Ax where the wave with frequency a/ 
occurs. We can see that these modes satisfy the relation k" + k' = 2k and u>" + a/ — 2u>. 


57.6 eV. 

In Fig. 0 it is evident that the frequency and wave number of the third mode satisfy 
u>" = 2u — u>' and k" = 2k — k', respectively. This mode can be visualized through the 
spectra of E y plotted in Fig. [6j where the mode at frequency uo" = 1.23a; manifest itself 
clearly. 

The evolution of E x spectrum in plasmas with higher electron temperature (123.9 eV) or 
weaker incident wave (E\ =0.25 MV/m) are also investigated. 

When temperature is higher, the location where the wave with frequency a/ = 0.77a; 
appears shifts to x ~ 680Ax, which is much deeper in the plasma. This means that the 
wave has a greater tendency to penetrate the plasma with higher temperature. At the same 
time, the reflectivity is lower. 

If we use a smaller source amplitude E \, the waves at frequencies a/ and to" generated 
during the mode conversion are much weaker than the wave with frequency a;. This is of 
course due to the nonlinear nature of the PDI. 
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FIG. 6: Frequency spectra of E y (x,t). The scenario is the same as Fig. [4] 


C. Dependence of reflectivity on wave amplitude 

In wave heating and current drive the energy conversion efficiency C is a key index. Given 
the energy reflectivity R, the conversion efficiency can be calculated according to C — 1 — R. 

During the nonlinear X-B mode conversion, the modes with different frequencies coexist. 

n 

The reflectivity defined for linear case, as in Ref. [14J, is no longer valid. In this work, we 
calculate R based on the Poynting vectors (see Appendix [B] for detail). 

We now consider the reflectivity of two incident electromagnetic waves with different 
amplitudes, i.e., Ei = 1 MV/m and Ei = 0.25 MV/m. The total number of time steps 
is 84000, the evolution of reflectivities with time are recorded in Fig. [TJ The reflectivity 
of the wave with the smaller amplitude (E\ = 0.25 MV/m) agrees well with the linear 
mode conversion theory at the beginning, but increases by 10% after sufficiently long time 
(t = 1200/u;). For stronger incident wave with amplitude E\ — \ MV/m, the reflectivity 
increases much faster. The reflectivity increases to almost 4 times of its original value after 
t = 1200/u. 

The simulation results show that the reflectivity is larger when the wave amplitude is 
higher. So in realistic wave heating scenario, incident waves with lower power density will 
have larger heating efficiencies. 
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FIG. 7: The evolution of reflectivities when waves with different amplitudes were injected 
into the non-uniform plasma. Here n r = 380 and T e = 57.6 eV. 

D. The energy deposition profile 


Since the location of wave energy deposition is crucial to rf heating, it is important to 
study the energy deposition profile. Generally speaking, the energy depositing inside the 
plasma core region will be helpful for heating and current driving. In this subsection, the 
energy deposition is obtained by comparing the kinetic energy distribution of particles at 
different times. In plasmas, the wave energy is stored in both electromagnetic fields and 
oscillating particles during each period of the wave. To calculate the energy deposition rate, 
we average the energies of particles over 4000 time steps and record the energy evolution in 
Fig. [8j 


According to Fig. 8(a) the energy deposition happens in a broad range in the plasma if 
the incident wave has a larger amplitude. Most wave energy is deposited near the location 
where the SHG of EBW appears (x ~ 608As). Different from the prediction of the linear 


theory, both Figs. 8(a) and 8(b) indicate that the wave energy penetrates to a deeper position 


instead of accumulating at x = xuhr exactly. These results can be used for the purpose of 
increasing the efficiency of auxiliary heating. 


The result in Fig. 8(b), which shows no evident electron heating, is consistent with the 
prediction of the linear theory. According to the linear theory, in the cold plasma approxi¬ 
mation the wave energy should accumulate at x = xuhr exactly. But in hot plasma, it is 
expected that the wave energy will be converted into the Bernstein mode and be carried into 
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(a) 


(b) 


FIG. 8: The difference of electron kinetic energy distribution at the beginning, middle, and 
end of the simulation. The amplitude of wave sources are 1.0 MV/m in subfigure (a) and 
0.25 MV/m in subfigure (b). To suppress the energy oscillation caused by the wave 
propagation, the kinetic energy are averaged over 4000 time steps. 


the plasma interior. For perpendicular propagation waves, the only linear damping of the 
Bernstein wave is caused by relativistic effects, which is weak for the electron temperatures 
in the simulation. So it is difficult for the wave to heat electrons. 


On the other hand, for waves with large amplitude in Fig. 8(a), obvious energy deposition 
occurs to the right of the UHR position, which is caused by the nonlinear effect. Simulation 
results show that a large amount of wave energy is accumulated at x ~ 608Ax, where the 


SHG of EBW appears. From Fig. 4(c) and Fig. 6(a) it can be observed that the waves 
with frequencies a/ = 0.77 oj and u" = 1.23a; are generated at x ~ 608Ax. Notice that 
to" = 1.23a; ~ 2a; ce (see Fig. EJ), and the wave with frequency to" will lead to electron 
heating through the second-harmonic resonance. The complete process of the nonlinear 
energy accumulation for large amplitude incident waves is provided as follows. 


1. The X wave converts to Bernstein wave with the same frequency a;. 

2. The Bernstein wave generates another Bernstein wave with frequency 2a; at x ~ 608Ax 
through nonlinear self-interaction. 

3. The new Bernstein wave with 2a; generates two other waves with frequencies a/ = 
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FIG. 9: The reflectivity of waves propagating in plasmas with different electron 
temperatures. Source amplitude E\ is 1.0 MV/m and n r = 380. 

0.77a; and uj" = 1.23a; ~ 2a; ce through PDI. 

4. The wave with uj" = 1.23a; ~ 2a; ce heats the electrons due to the second-harmonic 
resonance. 

Because the nonlinear effect strongly depends on the parameters of the plasma and incident 
wave, different nonlinear deposition processes will appear in different situations. We will 
investigate the large variety of nonlinear energy deposition phenomena in future studies. 

E. Reflectivity influenced by electron temperature 

To find out the dependence of the energy conversion efficiency on the background electron 
temperature, waves propagating in plasmas with different electron temperatures are investi¬ 
gated. The reflectivity evolution in plasmas with electron temperatures T e = 123.9 eV and 
T e = 57.6 eV is plotted in Fig. [9] We can observe that a long time (1200/a;) after the wave is 
injected, the energy reflectivity in the lower T e is about twice as large as the higher T e . This 
suggests that in wave heating experiments a higher electron temperature should be kept in 
order to maintain a larger mode conversion efficiency. 

To obtain more information of the relation between electron temperature and reflectivity, 
the case of T e ^0 eV is also tested. Two different wave amplitudes Ei = 0.066kV/m and 
Ei = 1.0 MV/m are studied. Fig. flOl shows that a short time (f = 156.484/a;) after the 
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FIG. 10: The spatial distribution E y (x) of wave injected from the source with small 
amplitude E\ = 0.066kV/m at t = 156.484/cu compared with linear theory. Here n r = 380 
and T e ~ 0. 



FIG. 11: The energy reflectivity evolution of waves with two different amplitudes for the 
case of T e « 0. 

wave is injected, E y (x) of the wave with small amplitude agrees well with the linear theory. 
On the other hand, Fig. |TT| implies that even for waves with small amplitude, the mode 
conversion structure is not stable and the reflectivity varies significantly after a long time 
(t fa 400/cj). Figure dU also suggests that X-B mode conversion in cold plasmas for wave 
with large amplitude (E\ =1.0 MV/m) is unstable, and the reflectivity increases quickly 
after t > 400 /u. 
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FIG. 12: Comparison of reflectivity with different density gradients. Dashed lines are the 
reflectivity obtained from linear theory. Here T e = 57.6 eV. 

F. Dependence of reflectivity on electron density gradient 

In this subsection, we investigate the dependence of reflectivity on the gradient of electron 
density. Thus we let n r in Eq. (jT8j) to be {380, 480, 600, 780}, T e = 57.6 eV, and E\ = 
1.0 MV/m. The reflectivities with different parameters are presented in Fig. [£2J We observe 
that under all the parameters, the reflectivities keep increasing during the mode conversion. 
The variation of reflectivity is very sensitive to the value of density gradient. Therefore, when 
applying the X-B rf heating in the magnetic confinement devices, it is preferable to maintain 
a large edge electron density gradient in order to maximize the conversion efficiency. 


IV. SUMMARY AND DISCUSSION 

The variational symplectic PIC scheme has been applied to simulate the nonlinear X-B 
mode conversion. The preservation of symplectic structure guarantees the accuracy of the 
long-term simulation results. Simulation results show that when the amplitude of incident 
wave is sufficiently large, the reflectivity will increase during the mode conversion, and 
modes with frequencies other than c o are exited. Even for waves with small amplitudes, the 
reflectivity changes significantly after a long time. In the meantime, we And that the energy 
deposits within a wide spatial range in the plasma. The simulation study also suggests that 
a higher electron temperature can improve the efficiency of energy conversion. Moreover, a 
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FIG. 13: Comparison of the reflectivity for two cases with different number of sampling 
particles per grid. It is evident that the physical results are not affected by the number of 
sampling particles. 


sensitive correlation between density gradient and reflectivity are observed. 

In PIC simulations, the numerical noise always exists due to the relatively less number 
of sampling particles. Large numerical noise may lead to inaccurate simulation results. The 
accuracy can be verified by checking the influence of the number of the sampling points on 
the simulation results. In the simulation, we have adjusted the number of sampling particles 
per grid from 4000 to 2000 and repeated the simulation with the same physical parameters 
to verify that the simulation results are not affected by the numerical noises. See Fig. [13] 

In this paper, the X-B mode conversion is discussed, and a series of interesting and useful 
nonlinear phenomena are investigated. The conservation property of the variational sym- 
plectic PIC plays a crucial role in achieving the new results during the long-term simulations. 
In future studies, the long-term simulations of mode conversion of other waves, such as lower 
hybrid wave (LHW) and ion cyclotron wave (ICW), will be carried out. We expect that, 
the variational symplectic PIC simulation is effective for these nonlinear rf processes and is 
able to suggest optimal parameters for rf heating and current drive schemes. 
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Appendix A: Linear Fluid Theory for the Reflectivity of Extraordinary Wave 


In fluid theory, the evolution of plasma is determined by the following equations, 


d . . 

—n e = -V • (n e u) , 


d 


m e ( — + u ■ V)u = e ( E + -u x B ) , 


1 f? _ _ ,, 47r , . 

E = V x B-(en e u) 

cot c 

~B = -VxE. 

c ot 


(Al) 

(A2) 

(A3) 

(A4) 


For the propagation of X-waves, we only need to consider the extraordinary component 
(Ei • B 0 = 0). For ID inhomogeneous, field variable can be expressed as 


/(r,t) = fo(x) + fi{x)e luJt 


After linearizing Eqs. ( 1Alj) - (lA4l) . we obtain a second-order ODE 

E"(o + g(o^(0 = o, 


(AS) 


(A6) 


where 


g(0 = i- 


Upe ~ ^p e ) 

(U 2 ((U 2 - (U 2 e - (U 2 e 


) ’ 


e = —, 

CCJ 

2 47rn e (a;)e 2 

^pe 

p rn e 


(AT) 

(AS) 

(A9) 


Here n e (x) is decided by Eq. (fTSjh n r = 380, and the external magnetic field is B 0 = 0.55 Te z . 
Note that this ODE may have a singular point £uhr where 


2 i Z Z A 

+ W ce “ W = 0 ' 


(A10) 


According to previous results 


of Refs, [loj] and fll] . 


the connection condition at £ = £ uhr is 


E ykuHR- o - e v\suhr+o ~ E v (€uhr) , (All) 

E y\(,UHR— o E y\iuHR+Q i7T (3Ey(£uHFt) i (A12) 

P = li m (£ “ £uhr) Q (£) • (A13) 

Z^ZUHR 

With these connection conditions, Eq. (1A6|) can be solved numerically. To find an appropri- 
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FIG. 14: The relation between amplitude reflectivity and cu, obtained by the linear fluid 
theory. Here n r = 380. 

ate <jj for the incident wave, we calculated the reflectivity for c o G [1.2 oo ce , 2.0cu ce ] and plotted 
in Fig.[[4j which shows that the lowest reflectivity is at uj ~ 1.37a ; ce . For u> G [1.2o; ce , 2.0cu ce ], 
the reflectivity will not be too high to make the mode conversion insignificant. For all our 
examples, u is set to 1.62cu C e, and the corresponding reflectivity is R ~ 0.15 when n r = 380. 

Appendix B: The Reflectivity Calculated by Averaged Poynting Vectors 

The new reflectivity is defined as follows. Consider the Poynting vector in the front of 
the source (x = 8Ax). The total energy flux S a ^ consists of incoming and outcoming parts, 
which can be expressed as 


(Bl) 



The incoming part is the energy flux of incoming wave 



(B2) 


where E[ is the amplitude of incoming wave. Note that E[ is not the amplitude of source 
E\. In our study, E[ = E\/2. The total energy flux is defined as 


S all = —E X B , 


(B3) 


and the outcoming part of energy flux is 



(B4) 


21 



Therefore, the energy reflectivity is 

R = ^ . (B5) 

^in 

All Poynting vectors here are time averaged in 207r/a;, which is significantly larger than the 
period of waves under investigation. 
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